16S Analysis Example
Rachel Rodgers
August, 23, 2022
knitr::opts_chunk$set(warning = FALSE,
message = FALSE)library("phyloseq")
library("ggpubr")
library("ggrepel")
library("data.table")
library("microbiome")
library("vegan")
library("gridExtra")
library("DESeq2")
library("plotly")
library("tidyverse")Read-In and Modify Data
A 16S analysis using phyloseq requires a phyloseq object (created from the DADA2 pipeline) and metadata, which describes each sample. Typically metadata is supplied via a mapping file which is created when you prepare your samples for sequencing.
First we need to read-in the sample metadata and the raw phyloseq data. In the modify-sample-data code chunk we can clean up the sample metadata as needed for the analysis. Then we’ll combine it with the phyloseq data. The phyloseq object is the basic object that we will use for most of the 16S analysis. See this phyloseq tutorial for more information.
sampleDataRaw <- read.delim("../documents/mappingFile_16S_Analysis_Workshop.txt",
check.names = FALSE)There are several columns in sampleDataRaw. Several of the columns aren’t required for this analysis. The columns that we really care about are the #SampleID and the Experiment column, which contains the variable of interest (which facility the sample came from). Let’s rename that column to something more reasonable, like Site.
sampleDataModified <- sampleDataRaw %>%
rename(Sample = `#SampleID`, Site = Experiment)As you go along in the analysis, you may realize that you need to modify the sample metadata in some way, such as adding or changing variables. You can return to the above code chunk and add code to modify as needed. This will keep all the sample data modifications in one location that is easily found, instead of little changes made here-and-there all throughout the script. Just remember to re-merge it with the phyloseq object, as shown in the next chunk.
physeqRaw <- readRDS("../data/physeqObjects/ps0.rdp_single.RDS")
physeqRaw # 508 taxa x 12 samples## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 508 taxa and 12 samples ]
## tax_table() Taxonomy Table: [ 508 taxa by 7 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 508 tips and 506 internal nodes ]# MAdd rownames to sampleDataModified before merging.
sampleDataFinal <- sampleDataModified
row.names(sampleDataFinal) <- sampleDataFinal$Sample
physeqMerged <- merge_phyloseq(sample_data(sampleDataFinal), physeqRaw)
physeqMerged## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 508 taxa and 12 samples ]
## sample_data() Sample Data: [ 12 samples by 10 sample variables ]
## tax_table() Taxonomy Table: [ 508 taxa by 7 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 508 tips and 506 internal nodes ]Read Count Distributions
At the beginning of the analysis, we check for sample outliers (samples with underlying data that do not conform to experimental or biological expectations) in order to minimize technical variance. There are several ways to detect these outliers. Here we will look for samples with unexpected numbers of read counts, then check how the samples cluster to see if samples are behaving unexpectedly.
Reads per Sample
Next we’ll look at how many reads are in each sample, and check if any samples have an unexpected number of reads. Some very low read samples may be low due to technical problems, as opposed to being due to true biological variation. Identifying (and removing) these types of samples helps us identify real variations in the microbiome due to biology and not some other confounder.
readsPerSample <- data.frame("reads_per_sample" = sample_sums(physeqMerged))
summary(readsPerSample)## reads_per_sample
## Min. : 7229
## 1st Qu.:12780
## Median :14373
## Mean :14657
## 3rd Qu.:17504
## Max. :21593# Merge readsPerSample back with the sample data for plotting.
sampleData <- merge(sample_data(physeqMerged), readsPerSample,
by = "row.names", all = TRUE)
sampleData <- column_to_rownames(sampleData, var = "Row.names")
# Let's sort the sampleData data frame from highest to lowest count sample.
sampleData <- sampleData %>%
rownames_to_column() %>%
dplyr::arrange(desc(reads_per_sample)) %>%
column_to_rownames()
sampleOrder <- as.character(sampleData$Sample)
# Now apply sampleOrder to the Sample column so they sort correctly in the plot
sampleData$Sample <- factor(sampleData$Sample, levels = sampleOrder)
# Plot the reads per sample bar plot
readsPerSampleBar <- ggplot(sampleData,
aes(x = Sample, y = reads_per_sample)) +
geom_bar(stat = "identity") +
xlab("Sample") + ylab("Number of Reads") +
ggtitle("Reads per Sample Bar Plot") +
#scale_y_log10() +
theme(axis.text.x = element_text(angle = 45),
plot.title = element_text(hjust = 0.5))
readsPerSampleBar
# Plot histogram showing distirubtion of read counts
readsPerSampleHist <- ggplot(sampleData, aes(reads_per_sample)) +
geom_histogram(color = "black", binwidth = 1000) +
xlab("Number of Reads") + ylab("Number of Samples") +
ggtitle("Reads per Sample") +
theme(plot.title = element_text(hjust = 0.5))
readsPerSampleHist
We can see here that there is not a big drop in read counts and no samples with unusually low read counts that may need to be removed. Also note that read counts have a somewhat normal (but not perfect) distribution.
Reads per Group
Since we will be checking for differences between the Site variable (Facility_1 vs. Facility_2) we should check for differences in the average read number. Significant differences here could have an impact on significant differences we find in ecological indices that we calculate, such as richness, so it’s good to know ahead of time.
readsPerGroup <- ggplot(sampleData, aes(x = Site, y = reads_per_sample)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(width = 0.2) +
xlab("Facility") + ylab("Reads per Sample") +
ggtitle("Average Reads per Facility") +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5)) +
stat_compare_means(method = "wilcox.test", label.x.npc = 0.5)
readsPerGroup
Statistically there is no difference in the average number of reads per sample by facility.
MDS Outliers
Creating an MDS plot of the Bray-Curtis distances between samples can give us some idea if samples are not clustering with other samples as may be expected. This could be indicative of some technical issue with the sample.
set.seed(32688024) # set the seed for reproducibility of "random" computations
brayOrd <- ordinate(physeqMerged, method = "MDS", distance = "bray")
mdsPlot <- plot_ordination(physeqMerged, brayOrd, color = "Site") +
theme_bw() +
geom_point(size = 3.5) +
ggtitle("MDS of Bray Distances") +
theme(plot.title = element_text(hjust = 0.5)) +
geom_text_repel(aes(label = Sample))
mdsPlot
While stool samples separate nicely by facility, there is also some sort of separation of samples from Facility_1 into two groups of 3. By looking at the bar plot above of reads per sample, it does not seem that samples are separating by read counts (low and high), so some other variable(s) may be at play. At this point there isn’t justification for removing any samples. We will explore distance metrics further in the analysis when we get to beta-diversity.
Taxa Filtering
Non-bacterial taxa need to be removed from the data set. We also want to check for low abunance/low prevalence taxa which may not contribute to the overall community evaluation or differential abundance testing.
# Note that this step will remove any taxa designated as NA at the Kingdom or
# Phylum rank.
# originally 508 taxa in 12 samples (physeqMerged)
physeqBacteria <- physeqMerged %>%
subset_taxa(Kingdom == "Bacteria" & Phylum != "Cyanobacteria/Chloroplast")
physeqBacteria # 498 taxa in 12 samples## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 498 taxa and 12 samples ]
## sample_data() Sample Data: [ 12 samples by 10 sample variables ]
## tax_table() Taxonomy Table: [ 498 taxa by 7 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 498 tips and 496 internal nodes ]Analysis can benefit from the removal of low prevalence and low abundance taxa. Prevalence is defined as the number of samples in which a taxon appears at least once. Here I will explore the relationship of prevalence and total read count for each feature at the phylum level.
# Calculate prevalence of taxa at a given taxonomic rank for low prevalence taxon filtering.
# Create a named vector where each element name is an OTU sequeunce, and each value
# is the number of sequences in which that OTU is present (max value is the total
# number of sequences).
prevalence_vector <- apply(X = otu_table(physeqBacteria),
MARGIN = ifelse(taxa_are_rows(physeqBacteria), yes = 1, no =2),
FUN = function(x) {sum(x > 0)})
# Generate a prevalence dataframe that also adds a TotalAbundance column (the total
# number of reads for that OTU across all samples) and the taxonomy information
# for each OTU.
prevalence_df <- data.frame(Prevalence = prevalence_vector,
TotalAbundance = taxa_sums(physeqBacteria),
tax_table(physeqBacteria))
phylaPrevalencePlot <- ggplot(prevalence_df,
aes(x = TotalAbundance,
y = Prevalence/nsamples(physeqBacteria), # to get a fraction
color = Family)) +
geom_hline(yintercept = 1/12, alpha = 0.5, linetype = 2) +
geom_point(size = 3, alpha = 0.7) +
scale_x_log10() +
xlab("Total Abundance") +
ylab("Prevalence [Frac. Samples]") +
facet_wrap(~ Phylum) +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5)) +
ggtitle("Phylum Prevalence in All Samples",
subtitle = "Colored by Familiy")
phylaPrevalencePlot
We see that Chlamydiae is has both low prevalence and low abaundance in this data set and may be uninformative. There are similar taxa within the Bacteroidetes, Firmicutes, Proteobacteria and Tenericutes. It is good to keep this in mind when we start analyzing differentially abundant taxa. If many taxa are returned through differential abundance testing, we may want to remove some of these low prevalent taxa to help hone in on which bacteria are truly meaningful in describing differences between our variable of interest.
Community Composition Plotting
Community composition plots can be used to examine overall taxon representation across samples. Here we’ll compare the taxa found in the samples from the different facilities, focusing at the Phylum level and without any abundance filtering.
abundanceDF <- physeqBacteria %>%
tax_glom(taxrank = "Phylum") %>%
transform_sample_counts(function(x) {x/sum(x)}) %>%
psmelt()
# If you wanted to add an abundance filter to, for example, remove any species
# that account for < 5% of the total reads in a sample, you could add after
# psmelt:
# %>% filter(Abundance > 0.05)
compositionPlot <- ggplot(abundanceDF, aes(x = Sample, y = Abundance, fill = Phylum)) +
geom_bar(stat = "identity", width = 1, color = "grey14") +
facet_wrap(~ Site, scales = "free") +
theme(axis.text.x = element_blank(),
axis.title.x = element_blank())
compositionPlot
You can make community composition plots at various taxonomic ranks and with different levels of abundance filtering. In a typical analysis, I would create a loop to make a composition plot at every level, or use an RShiny widget.
Alpha Diversity
Alpha diversity is a standard tool used to calculate ecological indices which give the number of taxa present in a study and the relationships between relative abundance and how evenly taxa are distributed. These indices provide useful summary information about the community structure of a given sample.
We are primarily interested in richness (the number of taxa per sample) and Shannon diversity (a description of the richness and evenness of a sample’s community). An evenness of 0 indicates that a community is dominated by one or a few taxa, but an evenness of 1 means that species are evenly distributed. Diversity increases as richness and evenness increase.
Here we’ll calculate richness (“observed”), Pielou’s evenness (“evenness_pielou”) and Shannon diversity (“diversity_shannon”), and add this information back to the sample data data frame for plotting.
alphaDiv <- microbiome::alpha(physeqBacteria,
index = c("observed", "evenness_pielou",
"diversity_shannon"))
sampleData <- merge(sampleData, alphaDiv, by = "row.names", all = TRUE)
sampleData <- column_to_rownames(sampleData, var = "Row.names")
# Make a boxplot of samples' richness and diversity, grouped by Site
richnessPlot <- ggplot(sampleData,
aes(x = Site, y = observed)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(width = 0.2) +
ylab("Richness") +
xlab("Facility") +
ggtitle("Richness") +
theme_pubr() +
theme(plot.title = element_text(hjust = 0.5)) +
stat_compare_means(method = "t.test", label.x.npc = 0.5)
evennessPlot <- ggplot(sampleData,
aes(x = Site, y = evenness_pielou)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(width = 0.2) +
ylab("Pielou's Evenness") +
xlab("Facility") +
ggtitle("Pielou's Evenness") +
theme_pubr() +
theme(plot.title = element_text(hjust = 0.5)) +
stat_compare_means(method = "t.test", label.x.npc = 0.5)
diversityPlot <- ggplot(sampleData,
aes(x = Site, y = diversity_shannon)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(width = 0.2) +
ylab("Shannon Diversity") +
xlab("Facility") +
ggtitle("Shannon Diversity") +
theme_pubr() +
theme(plot.title = element_text(hjust = 0.5)) +
stat_compare_means(method = "t.test", label.x.npc = 0.5)
grid.arrange(richnessPlot, evennessPlot, diversityPlot, ncol = 3)
Here we can see that while richness doesn’t vary by site, Shannon diversity is significantly lower for samples from Facility 2. This is likely being driven by the difference in evenness, which on its own is only nearly significant.
Beta Diversity
Distance Measures
While alpha diversity was concerned with describing the community within a sample, beta diversity compares the communities between samples. Beta diversity is quantified through “association” coefficients - similarity or distance measures. It describes how similar/dissimilar samples’ communities are to one another. We view these distances with ordination plots such as PCoA.
A common distance measure is the UniFrac distance, which is a measure that considers the phylogenetic relationships between bacterial taxa. A community of more closely-related taxa is less diverse than a community of distantly-related taxa. The weighted UniFrac measure adds on to this by accounting for the relative abundance of each taxa. Again, the choice of metric requires understanding of the underlying data, experiment and justification.
Here we will just focus on the weighted UniFrac distances.
To determine whether clustering is significant by Site, we will run an ADONIS test.
set.seed(46395617)
ord_wUniFrac <- ordinate(physeqBacteria,
method = "PCoA",
distance = "wunifrac")
pcoa_wUniFrac <- plot_ordination(physeqBacteria,
ord_wUniFrac,
color = "Site") +
theme_bw() +
geom_point(size = 4) +
ggtitle("PCoA of Weighted UniFrac Distance",
subtitle = "Colored by Facility") +
stat_ellipse(type = "norm") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
pcoa_wUniFrac
set.seed(46395617)
bdist <- phyloseq::distance(physeqBacteria, "wunifrac")
col <- as(sample_data(physeqBacteria), "data.frame")[, "Site"]
adonis.bdist <- adonis(bdist ~ col)
adonis.bdist$aov.tab## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## col 1 0.28476 0.28476 18.984 0.65498 0.002 **
## Residuals 10 0.15000 0.01500 0.34502
## Total 11 0.43476 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1According to ADONIS, the grouping of samples by Site is statistically significant. The R2 value indicates that approximately 66% of the variation in distances is explained by this grouping (whether a sample came from facility 1 or facility 2) and is significant when explaining the clustering of samples.
Differentially Abundant Taxa (Biomarker Analysis)
In this step, we are looking for differentially abundant taxa that may diagnose whether a sample originated from facility 1 or facility 2. A differentially abundant taxon is a species whose abundance is significantly different between categories. While there are many ways to test for these taxa, here we are using DESeq2.
dds <- phyloseq_to_deseq2(physeqBacteria, design = ~ Site)
# Run analysis
ddsAnalysis <- DESeq(dds, test = "Wald", fitType = "local", betaPrior = FALSE)
# Extract and format results
ddsResults <- results(ddsAnalysis,
contrast = c("Site", "Facility_1", "Facility_2"))
mcols(ddsResults)## DataFrame with 6 rows and 2 columns
## type description
## <character> <character>
## baseMean intermediate mean of normalized c..
## log2FoldChange results log2 fold change (ML..
## lfcSE results standard error: Site..
## stat results Wald statistic: Site..
## pvalue results Wald test p-value: S..
## padj results BH adjusted p-values# This contrast (and the mcols() results) tells us that results should be interpreted
# as facility 1 compared to facility 2. So, a positive log2FC means the taxon
# is more abundant in facility 1 and less abunant in facility 2. The opposite
# is true for a negative log2FC.
# From the DESeq results generated by GenerateDESeqResults, create a
# results data table that includes the taxonomy information and a column
# indicating whether results for each taxon are significant.
# Extract taxonomy table:
taxTable <- data.table(data.frame(as(tax_table(physeqBacteria), "matrix")),
keep.rownames = TRUE)
setnames(taxTable, "rn", "OTU")
setkeyv(taxTable, "OTU")
# Extract DESeq results as a data frame:
resDT <- data.table(as(ddsResults, "data.frame"),
keep.rownames = TRUE)
setnames(resDT, "rn", "OTU")
setkeyv(resDT, "OTU")
# Combine taxonomy information with the results table:
resDT <- taxTable[resDT]
resDT <- resDT %>%
filter(padj != "NA") %>%
mutate(Significant = padj < 0.05)
# Create volcano plot object
volcano <- ggplot(resDT,
aes(x = log2FoldChange, y = -log10(padj),
label1 = Family, # These labels are for plotly
label2 = Genus,
label3 = Species)) +
geom_point(data = subset(resDT,
resDT$Significant == FALSE), # non sig. taxa are grey
color = "grey") +
geom_point(data = subset(resDT,
resDT$Significant == TRUE),
aes(color = Phylum, size = baseMean)) + # sig. taxa colored by Phylum
geom_vline(xintercept = 0, lty = 2) +
geom_hline(yintercept = -log10(0.05)) +
ggtitle("Differntially Abundant Taxa by Facility") +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(angle = -90, hjust = 0, vjust = 0.5),
axis.title = element_text(size = 12),
axis.text = element_text(size = 12),
legend.text = element_text(size = 12))# save your script before interacting with the plotly plot. It has a tendency to
# crash the RStudio session.
# If you have problems, you can view the static plot directly by typing volcano.
ggplotly(volcano, tooltip = c("Phylum", "Genus", "Species", "log2FoldChange", "baseMean"))On initial glance, it appears that Bacteroidetes as a whole are more diagnostic of samples originating from facility 2, while Firmicutes are seen more in samples from facility 1. If you check back with the community composition plot, you can see that this pattern looks mostly true - Bacteroidetes appear to make up a larger portion of the overall sample from facility 2 compared to a sample from facility 1. And facility 1 samples appear to have higher levels of Firmicutes than facility 2 samples.
Program Info & System Information
This info is good for reproducibility and debugging.
Sys.Date()## [1] "2022-08-23"getwd()## [1] "C:/Users/rache/Box Sync/Computational_Info/16S_Analyses/16S_Analysis_Workshop/scripts"sessionInfo()## R version 4.0.5 (2021-03-31)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 16299)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] forcats_0.5.1 stringr_1.4.0
## [3] dplyr_1.0.9 purrr_0.3.4
## [5] readr_2.1.2 tidyr_1.2.0
## [7] tibble_3.1.8 tidyverse_1.3.2
## [9] plotly_4.10.0 DESeq2_1.30.1
## [11] SummarizedExperiment_1.20.0 Biobase_2.50.0
## [13] MatrixGenerics_1.2.1 matrixStats_0.62.0
## [15] GenomicRanges_1.42.0 GenomeInfoDb_1.26.7
## [17] IRanges_2.24.1 S4Vectors_0.28.1
## [19] BiocGenerics_0.36.1 gridExtra_2.3
## [21] vegan_2.6-2 lattice_0.20-41
## [23] permute_0.9-7 microbiome_1.12.0
## [25] data.table_1.14.2 ggrepel_0.9.1
## [27] ggpubr_0.4.0 ggplot2_3.3.6
## [29] phyloseq_1.34.0
##
## loaded via a namespace (and not attached):
## [1] readxl_1.4.1 backports_1.4.1 plyr_1.8.7
## [4] igraph_1.3.0 lazyeval_0.2.2 splines_4.0.5
## [7] crosstalk_1.2.0 BiocParallel_1.24.1 digest_0.6.29
## [10] foreach_1.5.2 htmltools_0.5.3 fansi_1.0.3
## [13] magrittr_2.0.3 memoise_2.0.1 googlesheets4_1.0.1
## [16] cluster_2.1.1 tzdb_0.3.0 Biostrings_2.58.0
## [19] annotate_1.68.0 modelr_0.1.9 colorspace_2.0-3
## [22] blob_1.2.3 rvest_1.0.3 haven_2.5.0
## [25] xfun_0.32 crayon_1.5.1 RCurl_1.98-1.6
## [28] jsonlite_1.8.0 genefilter_1.72.1 survival_3.2-10
## [31] iterators_1.0.14 ape_5.6-2 glue_1.6.2
## [34] gtable_0.3.0 gargle_1.2.0 zlibbioc_1.36.0
## [37] XVector_0.30.0 DelayedArray_0.16.3 car_3.1-0
## [40] Rhdf5lib_1.12.1 abind_1.4-5 scales_1.2.1
## [43] DBI_1.1.3 rstatix_0.7.0 Rcpp_1.0.9
## [46] viridisLite_0.4.1 xtable_1.8-4 bit_4.0.4
## [49] htmlwidgets_1.5.4 httr_1.4.4 RColorBrewer_1.1-3
## [52] ellipsis_0.3.2 farver_2.1.1 pkgconfig_2.0.3
## [55] XML_3.99-0.9 sass_0.4.2 dbplyr_2.2.1
## [58] locfit_1.5-8 utf8_1.2.2 tidyselect_1.1.2
## [61] labeling_0.4.2 rlang_1.0.4 reshape2_1.4.4
## [64] AnnotationDbi_1.52.0 munsell_0.5.0 cellranger_1.1.0
## [67] tools_4.0.5 cachem_1.0.6 cli_3.3.0
## [70] generics_0.1.3 RSQLite_2.2.12 ade4_1.7-19
## [73] broom_1.0.0 evaluate_0.16 biomformat_1.18.0
## [76] fastmap_1.1.0 yaml_2.3.5 knitr_1.39
## [79] bit64_4.0.5 fs_1.5.2 nlme_3.1-152
## [82] xml2_1.3.3 compiler_4.0.5 rstudioapi_0.13
## [85] ggsignif_0.6.3 reprex_2.0.2 geneplotter_1.68.0
## [88] bslib_0.4.0 stringi_1.7.8 highr_0.9
## [91] Matrix_1.3-2 multtest_2.46.0 vctrs_0.4.1
## [94] pillar_1.8.1 lifecycle_1.0.1 rhdf5filters_1.2.1
## [97] jquerylib_0.1.4 bitops_1.0-7 R6_2.5.1
## [100] codetools_0.2-18 MASS_7.3-53.1 assertthat_0.2.1
## [103] rhdf5_2.34.0 withr_2.5.0 GenomeInfoDbData_1.2.4
## [106] mgcv_1.8-34 hms_1.1.2 grid_4.0.5
## [109] rmarkdown_2.15 carData_3.0-5 googledrive_2.0.0
## [112] Rtsne_0.16 lubridate_1.8.0